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ABSTRACT 

A variety of molecular dynamics simulations of 
energetic materials is presented, demonstrating the ability 
to predict structural and thermodynamic properties of 
these materials. The studies are also used to explore, at 
an atomic level, dynamic processes that might influence 
conversion of the material to products. These studies are 
are presented to illustrate how information generated 
through molecular dynamics simulations can be used in 
the design, development and testing of energetic 
materials. 

1. INTRODUCTION 

The chemically interesting features of energetic 
materials have been advantageously employed in a wide 
variety of industrial and military applications, but often 
these utilizations have not been fully optimized. This was 
mainly due to the inability to identify and understand the 
individual fundamental physical and chemical steps that 
control the conversion of the material to its final products. 
The conversion of the material is usually not the result of 
a single-step reaction, or even a set of a few simple 
consecutive chemical reactions. Rather, it is an extremely 
complex process in which numerous chemical and 
physical events occur in a concerted and synergistic 
fashion, and whose reaction mechanisms are strongly 
dependent on a wide variety of factors. Direct 
measurements of mechanistic details that would provide a 
fundamental description of the conversion process are 
lacking due to substantial experimental obstacles. These 
difficulties have required the development of innovative 
theoretical methods and models designed to probe details 
of the various phenomena associated with the conversion 
of energetic materials to products. Toward this end, we 
have expended considerable effort in developing and 
critically assessing a realistic generalized model for use in 


molecular simulations of dynamic processes of condensed 
phase explosives. 

Our development of the model follows an 
evolutionary approach. We first start with a simple 
description of interatomic interactions between molecules, 
and apply the model in condensed phase molecular 
simulations in which each molecule is treated as a rigid 
entity. Several studies were performed to explore the 
ability of the model to reproduce structural and 
thermodynamic information, and to determine limits of 
this model and the rigid-molecule assumption when 
applied to various classes of CHNO explosives. Results 
indicate that within the low-pressure, low temperature 
regime, such an approximation is adequate for predicting 
structural and thermodynamic information. 

The next stage in our model development is to 
incorporate flexibility (not reaction) of the molecules into 
the model. Several simulations of increasing complexity 
have been performed to assess this extended model, and 
are summarized herein. These include investigations of 
nitromethane over the entire temperature ranges of both 
solid and liquid phases, and over large pressure ranges (0- 
14 GPa) in both phases, melting, and vibrational energy 
relaxation (VER) in liquid nitromethane after excitation 
of the C-H stretching vibrations. 

Our ultimate goal is to simulate conversion of the 
material to products at conditions of extreme temperatures 
and pressures; our next step in the model development 
will incorporate chemical reactivity. However, it is 
imperative that the chemical and physical events 
immediately preceding the conversion be accurately 
depicted. In this paper, we report numerous studies in 
which we have assessed such depictions, and describe 
lessons learned in the development of this generalized 
model of CHNO explosives. 
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2. DETAILS OF THE CALCULATIONS 
2.1. Molecular Simulation Methods 

The classical molecular simulation methods of 
molecular dynamics (MD) and molecular packing (MP) 
have been used to study the static and dynamic properties 
of energetic materials. These methods are limited by the 
classical approximation and an accurate description of the 
potential energy surface (PES) for the system. An MD 
simulation involves integrating equations of motion to 
generate temporal profiles of atomic positions and 
velocities, thus providing a dynamic description of the 
system. Also, thermodynamic information can be 
obtained by averaging properties evaluated at each 
integration step over the duration of the simulation. 
Molecular packing is an atomistic simulation method that 
is used to investigate structural features and properties 
very near local minima on the PES, and involves 
minimization of crystalline lattice energies by varying 
crystallographic parameters. MP cannot produce dynamic 
information; rather it provides information about 
equilibrium structures. 


The simple model, which does not have a description 
of intramolecular motions such as bond stretches, angle 
bends or torsional motions, can be used only in 
simulations in which the molecule is rigid. This model is 
suitable for calculations of simple thermodynamic and 
structural properties for low temperature and pressure 
regimes, where deformation of the molecule is not 
important. Assuming the rigid-molecule approximation 
significantly reduces computational expenses in 
simulations due to the elimination of costly terms in the 
interaction potential that describe molecular flexibility. 
However, regimes of extreme temperatures and pressures 
are of great importance in the study of energetic materials, 
and inclusion of molecular flexibility is required in order 
to study processes within these regimes. Thus, a 
subsequent important step was the extension of the SRT 
intermolecular potential to include the full intramolecular 
(non-reactive) interactions. This has been done for the 
case of the nitroalkane explosive, nitromethane. The 
modifications were simple additions of intramolecular 
terms to describe stretching, bending and torsional 
motions. These terms were parameterized using quantum 
mechanical information of the isolated molecule. 


2.2. Potential Energy Surface (PES) 


3. RESULTS 


Our initial model, hereafter denoted as the SRT 
model after the original authors (Sorescu et al., 1997), 
assumed that the potential energy used in this study for a 
system of N molecules can be described as the sum of 
intermolecular interaction terms: 

. N N 

Y Total =iyy Yintermolecular ^ ^ 


2 «=i j =i 


The intermolecular potential i consists of the 
superposition of a pairwise sum of Buckingham (6-exp) 
(repulsion and dispersion) and Coulombic (C) potentials 
of the form: 


v ap ( r > A ap exp (~B a/j r) - C aj3 tr 6 (2) 

and 
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where r is the interatomic distance between atoms a and 
P, q a and qp are the electrostatic charges on the atoms, and 
e 0 is the dielectric permittivity constant of free space. 


The set of partial charges used in these calculations 
were determined through fitting these to the quantum- 
mechanically derived electrostatic interaction potential for 
an isolated molecule whose atoms are arranged in the 
experimental crystallographic arrangement. The 
remaining exponential-six parameters were adjusted to 
reproduce the experimental structure of the RDX crystal 
at ambient conditions (Sorescu et al., 1997). 


3.1 Rigid Molecule Simulations 

In addition to reproducing the crystal structure of 
RDX at ambient conditions using both MP and MD 
techniques (Sorescu et al., 1997), we found that this 
interaction potential could also describe the geometric 
parameters and lattice energies of different polymorphic 
phases of two other nitramine crystals: the polycyclic 
nitramine HNIW (CL-20) (Sorescu et al. 1998a) and the 
monocyclic nitramine HMX (Sorescu et al. 1998b). 
Further investigations exploring the limits of 
transferability of this interaction potential to other 
energetic molecular crystals were undertaken by 
performing molecular packing calculations for 30 
nitramine crystals (Sorescu et al. 1998c) and 51 non- 
nitramine CHNO crystals (Sorescu et al. 1999a). MP 
calculations using this interaction potential reproduced the 
crystal structures of all of these to within 5% of 
experiment. An extremely important result of the MP 
studies was the sensitivity of the results upon the selection 
of the partial charges used in the study. As discussed in 
Section 2.2, these charges have been determined from fits 
to ah initio electrostatic potentials calculated for the 
individual molecules whose atoms are arranged in the 
experimental configurations We have considered four 
different electrostatic models, with charges determined at 
a variety of levels of quantum mechanical theory: 
Hartree-Fock (HF), gradient-corrected non-local Density 
Functional Theory using the B3LYP density functional 
(B3LYP), second order second-order Moller-Plesset 
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perturbation theory (MP2) and charges calculated at the 
HF level uniformly scaled by a factor of 0.9. There is 
only a small influence (generally less than 1%) on the 
crystallographic parameters by the set of electrostatic 
charges used. However, the lattice energies of the 
crystals are significantly influenced by the electrostatic 
model. In particular, the best agreement with the 
experimental lattice energies has been obtained for the 
MP2 charges. The lattice energies calculated using the 
B3LYP charges overestimate the MP2 energies by about 
2.6% while the HF charges overestimate the MP2 
energies by 13.6%. The procedure of uniformly scaling 
the HF charges by the 0.9 factor decreases the above 
differences to about 6.2%. Lattice energies calculated for 
the HMX, CL-20, TNT and PETN systems all support 
polymorphic stability rankings determined experimentally 
(Sorescu et al., 1998a, 1998b, 1998c, 1999a). 

The next step in model development was to explore 
the validity of the rigid-molecule assumption in 
simulations of energetic materials under hydrostatic 
compression. Therefore, we analyzed the dynamics of the 
energetic crystals RDX, HMX, HNIW and PETN under 
hydrostatic compression conditions using isothermal- 
isobaric (NPT) MD simulations and this simple 
intermolecular potential (Sorescu et al., 1999b). In that 
study, predicted lattice parameters for the RDX, HMX 
and HNIW crystals were found in good agreement with 
experimental values over the entire range of pressures 
investigated experimentally. For the PETN crystal, the 
calculated crystallographic parameters were in acceptable 
agreement with experimental data only for pressures of a 
few GPa. For higher pressures, the disagreements of 
predictions with experiment demonstrated the inadequacy 
of the rigid-body approximation when used in simulations 
of floppy molecules such as PETN. 

This effect is illustrated in Figure 1, in which the 
calculated volumes of the RDX and PETN crystals are 
plotted as functions of pressure, and compared with 
experiment. While the MP predictions, which do not 
include thermal effects, are in some disagreement with 
experiment, the NPT-MD simulations of RDX are in very 
good agreement with experiment over the range of 
pressures explored. However, the failure of the model to 
predict crystal densities of PETN under high compression 
is quite apparent. Therefore, the SRT results suggest that 
at moderate temperatures and pressures, simulations using 
the rigid-molecule approximation will provide reasonably 
accurate results at a significantly reduced computational 
cost compared to those that use more complex flexible 
interaction potentials. 

3.2 Structural studies of flexible nitromethane in the 

solid and liquid phase 


Sorescu, Rice, and Thompson extended the original 
SRT model by adding an intramolecular part which 
consists of a superposition of bond stretching, bond 
bending and torsional angle terms (Sorescu et al., 2000). 
In particular, Morse potentials have been used to 
represent bond stretches while harmonic and cosine type 
of potentials have been used to simulate the bending and 
torsional motions. These terms were parameterized based 
on the geometric and vibrational frequencies and the 
corresponding eigenvectors data obtained from ab initio 
molecular orbital calculations for the isolated molecule. 
Molecular packing calculations using the proposed 
potential produce an accurate prediction of the 
crystallographic parameters with deviations less than 
1.2% for the lattice edges. Moreover, NPT-MD 
simulations performed over the temperature range 4.2-228 
K and pressure range 0.3-7.0 GPa indicate that the 
crystallographic parameters are well reproduced for the 
entire range of temperatures and pressures simulated. 
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Figure 1. Comparison of the dependence of volume 
compression V/Vo on the external pressure obtained in 
MP and NPT-MD simulations for the a-RDX crystal 
(upper frame) and PETN (lower frame) with the 
experimental results. 


Excellent agreement was found for the calculated 
bulk modulus of nitromethane (6.78 GPa) with the 
experimental data (7.0 GPa). The corresponding 
predicted values as a function of either temperature or 
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pressure reproduce very well the experimental data 
obtained by neutron diffraction techniques (Trevino and 
Rymes, 1980) or by X-ray diffraction (Cromer et al., 
1985). Further, the data agree very well with the 
experimental findings in which the methyl group was 
found to be rotated by about 45° relative to the low 
temperature configuration. Figure 2 shows the time 
averaged distributions of the three Fl-C-N-0 dihedral 
angles in nitromethane, averaged over all molecules in the 
simulation cell, as functions of temperature and pressure. 
The distributions show that at 1 atm for all temperatures 
studied, the orientation of the methyl group oscillates 
about its equilibrium position determined at 4.2 K, latm. 
Flowever, the peaks of the distributions shift with 
increases in pressure at T=293 K. At low pressures, the 
distributions indicate that the orientation of the methyl 
groups are the same as that of the low-temperature, low- 
pressure crystal. There is a continuous shift of the peak 
positions with pressure such that between 0.3 GPa and 5.4 
GPa this shift amounts to about 41° while between 0.3 
GPa and 7.0 GP the corresponding variation is about 50°. 
Also, the corresponding activation energy for methyl 
rotation was found to be in the range of the reported 
experimental activation energies. 
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Figure 2. Distribution of the Flj-Ci-N 2-03 (i=5,6,7) 
dihedral angles for all nitromethane molecules in the 
simulation box as function of temperature (1 atm) (a) and 
pressure (T=293 K) (b). 

In a subsequent study, the transferability of the general 
intra- and inter-molecular potential developed for 
crystalline nitromethane to the liquid phase of 
nitromethane was explored by computing various 
physical properties of the liquid as functions of 
temperature and pressure (Sorescu et al., 2001). A large 
set of static and dynamic properties of liquid nitromethane 
have been considered in these tests including the heat of 


vaporization, the variation of density with temperature 
(over the range 255-374 K) and pressure (over the range 
0-14.2 GPa), the thermal expansion coefficient, the self¬ 
diffusion coefficients, the viscosity coefficient, the 
dielectric constant, the bulk modulus, and the variation of 
vibrational frequencies with pressure. The analyses 
performed using NPT-MD simulations show that the great 
majority of these structural, energetic and spectroscopic 
are well reproduced. The only exception is the dielectric 
permittivity, which was underestimated. This limitation 
was attributed to the lack of polarization effects in the 
intermolecular interactions. 

3.3 Melting of nitromethane 

Melting of nitromethane was also explored using the SRT 
model using two types of MD simulations. The first type 
of simulation is one in which the crystal is gradually 
heated until a parameter that monitors the degree of 
translational order in the crystal abruptly decreases. This 
change indicates that the system has transitioned from the 
crystalline to the liquid state, and the temperature at 
which this occurs is the “transition” temperature (see 
Figure 3). 

MD simulations using this method for several atomic 
crystals have shown that the transition temperature for a 
perfect crystal is substantially higher than the true 
thermodynamic melting temperature, but that the 
introduction of a critical number of voids lowers the 
transition temperature near the true thermodynamic 
melting point (Solca et al., 1997, 1998; Agrawal et al., 
2003b). Therefore, the key assumption in this method is 
that the true thermodynamic melting temperature for any 
model corresponds to the transition temperature for a 
crystal containing a critical concentration of voids. The 
calculated value of the melting temperature for the 
nitromethane model using this method is in good 
agreement with experiment; however, because this 
method of predicting the thermodynamic melting point is 
empirical, a second melting simulation was performed in 
which a coexistence of the liquid and solid phases were 
simulated to confirm the result. 

For this method, a simulation cell was constructed in 
which a block of liquid nitromethane was appended to a 
rectangular block of crystalline nitromethane, and the 
system equilibrated (using NPT-MD) to the desired 
temperature and pressure. Once equilibrated, NVE-MD 
simulations were performed and the temperature and 
behavior of the system monitored. If the temperature of 
the system is too high, the solid portion of the crystal will 
melt. If the temperature of the system is below the 
melting point, the liquid portion of the cell will solidify. 
The melting point is the temperature for which the liquid 
and solid maintain a coexistence. The results for the two 
methods were in near agreement, with the slight 
difference being attributed to hysteresis associated with 
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the direct heating process imposed in the void-nucleated 
melting simulation. 




Figure 3: Snapshot of nitromethane crystal resulting 
from molecular dynamics simulations before (top) and 
after (bottom) reaching the melting point. 



Figure 4: Pressure (in kbar) versus the melting 
temperature T (in K). The circles denote computed 
melting points, the line is a fit of the computed points to 
the Simon-Glatzel equation, the rectangles are 
experimental melting points reported by Jones and 
Giauque (Jones and Giauque, 1947) and the triangles 
denote the experimental values of Piermarini et al. 
(Piermarini et al., 1989). 

3.4 Energy transfer in liquid nitromethane 

Non-equilibrium molecular dynamics simulations 
were used to study vibrational energy relaxation (VER) in 
liquid nitromethane after excitation of the C-H stretching 
vibrations. This study was designed to explore the role 
of multiphonon up pumping in shock initiation of 
energetic materials, a theory based on mechanisms of 
energy flow in a shocked system. The multiphonon up 
pumping process begins with the heating of the phonons 
of a material upon the passage of a shock wave. The 
excess energy of the phonons flows into the cold 
molecular modes of the material through low-frequency 
modes that are strongly coupled to the phonon bath, and 
whose frequencies are near that of the maximum 
frequency of the phonon continuum (Chen et al., 1994). 


The values of the calculated melting temperature, T m , 
are found to be in good agreement with the experimental 
data at various values of pressure ranging from 1 atm to 
30 kbar (Figure 4). The computed values of the melting 
temperature satisfy the Simon-Glatzel equation: P(kbar) = 
aT m b + c, where a = 1.597xl0" 5 , b = 2.322, c = -6.74, and 
T m is in Kelvin. A comparison of computed T m with and 
without the presence of molecular vibrations reveals that 
T m is insensitive to the intramolecular interaction term of 
the potential energy function, but depends strongly on the 
intermolecular interactions, particularly the Coulombic 
term (i.e., the partial charges on atoms). 


As energy continues to flow from the overheated 
phonon bath into the doorway vibrations, higher-energy 
vibrational modes are subsequently excited through 
intramolecular vibrational energy transfer. This “up 
pumping” of vibrational states continues until sufficient 
energy is localized in the reaction coordinates of the 
molecules to allow reactions. Dlott and co-workers have 
performed extensive experimental investigations of 
energy transfer in condensed phase molecular systems to 
assess various aspects of this theory (Chen et al., 1994, 
1995; Hong et al., 1995; Deak et al. 1999; Dlott, 2001). 
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Dlott and co-workers have performed a series of 
spectroscopic studies of vibrational energy relaxation 
(VER) in liquid nitromethane (NM) that have 
demonstrated that the VER process in this system is fairly 
complex (Chen et al., 1994, 1995; Hong et al., 1995; 
Deak et al. 1999; Dlott, 2001). The most recent study 
(Deak et al., 1999) used anti-Stokes Raman spectroscopy 
to measure instantaneous vibrational populations of NM 
after infrared excitation in the C-H stretching region (near 
2970 cm" 1 ) for the neat liquid. Deak et al. (Deak et al., 
1999) report that this pulse excites the C-H stretching 
vibrations and the first overtones of the antisymmetric 
CH 3 bending and N0 2 stretching vibrations. Additionally, 
this IR-Raman method was applied to solutions of NM- 
CC1 4 to extract details of energy flow that could not be 
obtained from experiments on the neat liquid. CC1 4 is 
transparent at 2970 cm' 1 ; thus, it was used as a 
“molecular thermometer” to monitor the excitation of the 
bath upon vibrational cooling (VC) of NM. Monitoring 
the populations of CC1 4 vibrations after pump pulse 
excitation of the NM provided insight into the mechanism 
of VC for NM. The data show that the relaxation takes 
place in three steps. First, energy deposited in the C-H 
stretch (and in the first overtones of the antisymmetric 
NO 2 stretching and CH bending vibrations) is 
redistributed to all other vibrations within a few 
picoseconds. The NM-CCI 4 results do not reflect 
excitation of the CCI 4 during this time, indicating that the 
initial VER is intramolecular. Subsequently, the higher 
energy vibrations of NM (1560 and -1400 cm" 1 ) relax on 
the time scale of -15 ps, mainly by populating the lower 
energy vibrations (all transitions below -1400 cm" 1 ). 
Approximately one-third of the energy from the decay of 
these two transitions is dissipated to the bath. Finally, 
the lower energy vibrations excited in the first two stages 
relax by heating the bath. VC of NM occurs on the 50 to 
100 ps time scale, with the response of the CC1 4 
vibrational transitions after excitation of the NM rising on 
the same time scale. This indicates that the excitation of 
the bath molecules occurs mainly through indirect 
intermolecular vibrational energy transfer (IVET) 
processes. 

In an attempt to simulate the Deak et al. experiments 
(Deak et al., 1999) using classical MD, we utilized 
projection methods to follow the energy flow from 
excited molecular vibrational modes. In order to generate 
a fully-equilibrated liquid at the conditions of the 
experiment before excitation, an NPT MD simulation at 
294 K 1 atm was performed. Next, a microcanonical MD 
simulation of liquid NM was performed, in which a 
percentage of the molecules in the equilibrated liquid 
were selected for mode-specific excitation. We did not 
identify information in the experimental papers that 
would allow us to quantify the populations of the various 
vibrational states that are excited from the mid-IR pulse. 
Since we did not have quantitative information on 


vibrational populations, we arbitrarily chose to excite 120 
molecules (25% of the total number in the simulation cell) 
with the excitation energy equipartitioned among the 
three CH stretching vibrations. Each excited molecule 
was given a total of 2.075 kcal/mol in the form of kinetic 
energy equally partitioned among the three CH stretching 
modes, thus introducing an excess energy corresponding 
to a temperature rise of 12.4 K, a value that is close to the 
experimental temperature jump resulting from the pump 
pulses in neat NM (-10 K). 



Figure 4. Average kinetic energies for normal modes (a) 
v(CH 3 ); (b) v a (N0 2 ); (c) 5 a (CH 3 ); (d) 5 S (CH 3 ); (e) v s (N0 2 ) 
(f) p(CH 3 ); (g) v s (CN); (h) 5 S (N0 2 ); and (i) p(N0 2 ) of the 
excited nitromethane molecules. 

If the energy initially given to each molecule was 
immediately redistributed equally among its vibrational 
degrees of freedom, then the kinetic energy of vibration of 
each mode would increase from 0.292 to 0.341 kcal/mol. 
The information shown in Figure 4 indicates that such an 
immediate and uniform redistribution does not occur. 
Immediate increases in kinetic energy in all vibrational 
modes upon relaxation of the C-H stretches is evident. 
However, the differences in the curves show that energy 
transfer into the modes occurs at different rates, indicating 
that VC occurs in stages. The decay of the C-H stretching 
vibrational modes is exponential; however, the CH 3 
symmetric stretch [Fig. 4(a)] has a VER lifetime of 2.5 ps, 
whereas the CH 3 asymmetric stretching mode has a VER 
lifetime that is -3 times longer. The CH 3 asymmetric 
bends [Fig. 4(c)], N0 2 symmetric stretch [Fig. 4(e)] and 
CH 3 rocks [Fig. 4(f)] achieve their maximum level of 
excitation of -0.37-0.38 kcal/mol almost immediately. 
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Energy flow from these vibrations modes occurs 
exponentially, with decay constants ranging from 11.9 to 
16.2 ps. Four vibrational modes in the mid- to low- 
frequency range reach their maximum level of excitation 
at ~4 ps. The higher-frequency CF 13 symmetric bend [Fig. 
4(d)] and N0 2 asymmetric stretch [Fig. 4(b)] attain a 
maximum kinetic energy of 0.34 kcal/mol at 4 ps, and the 
lower-frequency N0 2 symmetric bend [Fig. 4(h)] and 
rocking modes [Fig. 4(i)] attain a maximum energy of ~ 
0.32 kcal/mol at 3.8 ps. Energy transfer from these 
modes can also be described as an exponential decay, but 
the lifetime of the low-frequency modes is much shorter 
than those of the higher frequency modes. The CN 
stretch mode [Fig. 4(g)] reaches its maximum excitation 
at 10 ps; energy flow from this point in time can be 
described as an exponential decay with lifetime ~20 ps. 
Although the kinetic energy for this mode appears to 
have reached the expected equilibrated value at 40 ps, the 
energy flow back into the mode, suggesting that it has not 
fully relaxed. 

The profdes of the vibrational modes of the unexcited 
NM molecules as functions of time indicate that some 
VER from the excited NM to the bath occurs at the 
beginning of the VC process, though the degree of 
excitation is not great and the energy appears to be 
uniformly distributed among the vibrational modes. The 
behavior of the curves supports an indirect mechanism of 
IVET, characterized by Deak et al. (Deak et al., 1999) as a 
two-step process in which VC of the vibrationally hot NM 
excites the phonons in the liquid, which subsequently 
excite vibrations of the solvent molecules by multiphonon 
up-pumping. 

Overall, the results are in qualitative agreement with 
experimental measurements of VER in liquid 
nitromethane after mid-IR excitation in the C-Fl stretching 
region The simulation results indicate that the excited 
C-FI stretching vibrations deposit energy predominantly 
into the remaining vibrations in the molecules. These 
vibrations relax at different rates, resulting in a multi¬ 
stage vibrational cooling process for nitromethane, in 
agreement with experimental results. The excitation of 
vibrations of the surrounding unexcited molecules occurs 
through indirect rather than direct intermolecular 
vibrational energy transfer processes, also in agreement 
with experiment. The main discrepancy between the 
experimental results and our results on the effect on the 
bath is the time scale on which the heating of the bath 
occurs. The experiments showed that while some energy 
buildup in the bath occurs on the 15 ps time scale, the 
majority of heating occurs on the 50-100 ps time scale 
corresponding to the final step of VC. The simulations 
show heating of the bath begins immediately, with full 
equilibration by 60 ps. 


3.5 Shock Hugoniot of Nitromethane 

The shock Flugoniot of a material provides a 
characterization of the behavior of a shock wave in that 
material, and is often used to assess performance of an 
explosive. The Hugoniot states are those that satisfy the 
Hugoniot relation (Zeldovich and Raizer, 1966) 

H=E-E 0 -'/2(R + P o )(V 0 -V)=0 (1) 



Figure 5. Shock Hugoniot of nitromethane. Calculated 
values from NPT-MD simulations using the SRT model 
are compared with experimental data (Mader, 2002; 
Delpuech and Mend, 1983) and results of thermochemical 
calculations using CHEETAH (Fried et al., 1998). The 
shock pressure is plotted versus the specific volume. 

where E and V are the specific internal energy and 
volume, respectively, and P is the pressure. The term 
specific refers to the quantities normalized to unit mass 
and the subscript “o” denotes the quantity in the 
quiescent, unshocked material. 

The equation of state of nitromethane described by 
the SRT potential was calculated using NPT-MD 
simulations; these results were used to calculate its shock 
Hugoniot. Results are shown in Fig. 5, along with 
experimental data and results calculated using the 
conventional thermochemical code CHEETAH and two 
different Equations of State. As evident in the figure, the 
Hugoniot curves predicted by CHEETAH using standard 
EOS are significantly different than that predicted by the 
SRT potential and experiment, and cannot adequately 
describe this system. 
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CONCLUSIONS 

Overall, the above results indicate that a wide variety 
of properties can be predicted by using molecular 
simulations and the SRT model. We expect that further 
extensions of the model to describe other systems and to 
allow breaking and formation of bonds will result in full 
atomistic description of the initiation and conversion of 
real energetic materials to final products. 
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